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A FINITE-ELEMENT ANALYSIS OF FATIGUE CRACK CLOSURE 
By J. C. Newman, Jr. 

NASA Langley Research Center 
Hampton, Virginia 

ABSTRACT 

Experiments have shown that fatigue cracks close at positive loads during 
constant-amplitude load cycling. The crack-closure phenomenon is caused by 
residual plastic deformations remaining in the wake of an advancing crack tip. 
The present paper is concerned with the application of a two-dimensional, non- 
linear, finite-element analysis for predicting crack-closure and crack-opening 
stresses during cyclic loading. 

A two-dimensional finite-element computer program, which accounts for both 
elastic-plastic material behavior and changing boundary conditions associated 
with crack extension and intermittent contact of the crack surfaces under cyclic 
loading, has been developed. An efficient technique to account for changing 
boundary conditions was also incorporated into the nonlinear analysis program. 
This program was subsequently used to study crack extension and crack closure 
under constant-amplitude and two-level block loading. The calculated crack- 
closure and crack-opening stresses were qualitatively consistent with experi- 
mental observations. 

INTRODUCTION 

Until recently, fatigue-crack propagation was assumed to be directly 
related to the linear elastic stress-intensity factor [1] during cyclic loading. 
Implicit in this concept were the assumptions that only the tensile portion of 
the load cycle was effective in growing the crack, and that cracks close pre- 
cisely at zero load. Elber [ 2-b ] has shown experimentally that fatigue cracks 
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close at positive loads during zero-tension const ant -amplitude load cycling. He 
has indicated that fatigue-crack closure may be a significant factor in causing 
the stress-interaction effects on crack growth (retardation or acceleration) 
under general cyclic loading. He has also postulated that the crack-closure 
phenomenon is caused by residual plastic deformations remaining in the wake of 
the advancing crack tip. 

Kobayashi, Chiu, and Beeuwkes [5], and Anderson [6] using the finite- 
element method analyzed the steadily growing crack under a single monotonically 
increasing load. Since these investigations did not consider cyclic loading, 
crack-closure effects were not accounted for in the analysis. 

Newman and Armen [7], also using the finite-element method, analyzed an 
extending crack under cyclic loading and included the effects of crack closure. 
However, theirs was a general-purpose program which included a large number of 
options concerning the types of elements and the material stress-strain behavior 
and, consequently, required a large computer storage capability. 

In the present paper, a more efficient two-dimensional, nonlinear, finite- 
element analysis [8] is presented, and used to investigate crack extension and 
crack closure in a center-cracked panel under cyclic loading. The panel material 
was assumed to be elastic -perfectly plastic. Both constant-amplitude and two- 
level block loading were investigated. The crack-closure stresses, crack- 
opening stresses, displacements, and residual stress distributions near the crack 
tip were determined as functions of applied stress. 

SYMBOLS 

a Half-length of crack, m 

Matrix relating total strains to nodal displacements, 

Material crack-growth constants 
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Incremental applied load vector, N 
Incremental load vector, N 

Differential volume of triangular element, vP 
Incremental total strain vector 
Incremental elastic stress vector, N/m 

2 

Incremental "initial" stress vector, N/m 

2 

Elasticity matrix relating stress to total strain, N/m 
Relaxation parameter 

Spring stiffness in x- and y-direction, respectively , N/m 
Diagonal matrix containing spring stiffnesses, N/m 
Elastic stiffness matrix, N/m 
Applied load vector, N 
"Effective" plastic load vector, N 

Stress ratio (ratio of minimum to maximum applied stress) 

2 

Applied gross stress, N/m 

2 

Maximum gross stress, N/m 
Minimum gross stress, N/m 

2 

Crack-opening stress, N/m 
Generalized nodal displacement vector, m 
Displacements in x- and y-direction, respectively, m 
Total plate width, m 
Cartesian coordinates 
Incremental crack growth, m 

2/2 

Effective stress-intensity factor range, N/m 
Incremental number of cycles 
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2 

Normal stress acting in y-direction, N/m 

2 

Tensile yield stress, N/m 


FINITE -ELEMENT ANALYSIS 

The elastic-plastic analysis of the center-cracked panel (Fig. 1) employed 
the finite-element method and the initial-stress concept as described in [9]. 

The finite-element model, Figure 2, was composed of two-dimensional constant- 
strain triangular elements (unit thickness). In the initial-stress approach, 
the load-displacement relations for a discretized structure are written to 
include the effects of initial stresses, which are required in order to satisfy 
the yield criterion (von Mises ) for an elastic-plastic material. These initial 
stresses produce effective plastic-load vectors which are applied to all elements 
which have become plastic and which maintain the permanent plastic deformation 
on those elements while the external loads are applied. The governing matrix 
equations for a discretized structure are reviewed only briefly here to demon- 
strate how the material nonlinearity is accounted for and what is required to 
treat crack extension and intermittent contact or separation of the crack 
surfaces. 

Solution Procedure for Elastic-Plastic Structures 

The application of the finite-element method to problems involving linearly 
elastic materials is straightforward because the material properties are constant 
and only one solution is required to obtain displacements for the elastic struc- 
ture. However, for elastic -plastic structures the coefficients in the stiffness 
matrix are functions of loading. Thus, the displacements are usually obtained 
by applying pmai 1 load increments to the structure and either updating the coef- 
ficients of the stiffness matrix or applying an "effective" plastic-load vector 
after each load increment. 
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In general, the matrix equation which governs the response of a discretized 
structure under loads which cause plastic deformation is 

[K e J (Ujj 1 = {p} 1 + (QJxli (1) 

where [K e ] is the elastic stiffness matrix, tu) is the generalized nodal dis- 
placement vector, (P } is the applied load vector, and {Q} is the "effective" 
plastic-load vector which accounts for elements in a plastic state. In the 
initial-Btress method, the solution to an elastic-plastic problem is obtained 
by applying a series of small load increments to the structure until the desired 
load is reached ((p) 1 = {p}*“' L + {dP}). The superscript i in Equation (l) 
denotes the current load increment and i - 1 denotes the preceding increment. 
After each load increment an iterative process is required to stabilize the 
plastic-load vector. The subscript I in Equation (l) denotes the current 
iteration and I - 1 denotes the preceding iteration. During the i increment 
a purely elastic problem is solved and the increments in total strain {de) and 
corresponding elastic stress {da e } are computed from the displacements for every 
element. Because of the material nonlinearity the stress increments are not, in 
general, correct. If the correct stress increment for the corresponding strain 
increment is {do), then a set of body forces or plastic-load vectors (dQ) caused 
by the "initial" stress (do 0 ) (*= {do g } - (do)) is required to maintain the stress 
components on the yield surface. The correct stress increment {da} is computed 
from the equations given in [9]* The plastic load increments are computed from 

CdQ} • Yi f CB3 T {da 0 } dV (2) 

m-l* / n 


5 



where M is the total number of elements, [B] is the strain-displacement rela- 
tionship, and the superscript T denotes the matrix transpose. The integration 
is taken over the volume of each element and the summation is over all elements 
in the structure. For elements which are in an elastic state, (da°) = 0. The 
total plastic-load vector is then computed as 

{<*}* = + (dQ) (3) 

At the second stage of computation the new force system (Q)^ is added to the 
applied load vector and a new set of displacements is obtained. Again, some of 
the stresses are likely to exceed the yield criterion and a new set of plastic- 
load increments is computed. The iteration process is repeated until the change 
in the plastic-load vector (Eq. (2)) is sufficiently small (0. 1 percent of the 
final values was chosen here). Usually, 5 to 15 iterations are required to 
stabilize the plastic-load vector. However, for conf igurations which have large 
strain gradients, more iterations are required. For the cracked plate considered 
here and the particular element mesh used, 10 to 30 iterations were required. 

In order to reduce the number of iterations, a relaxation technique was incor- 
porated into the nonlinear analysis program by using the equation 

(Q}j = + g{dft) (k) 

where g is the relaxation parameter. Because the displacements from the pre- 
ceding increment or iteration are used to compute the plastic-load increment, 
the plastic-load vector is underestimated. Thus, the relaxation parameter is 
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used to increase the plastic-load vector and, consequently, increase the rate of 
convergence. For the finite-element mesh used here, the displacements were 
found to converge roughly twice as fast using g = 2 than using g = 1 (usual 
value for the initial-stress method). 

Solution Procedure for Changing Boundary Conditions 

As previously mentioned, the finite- element analysis of an extending crack 
under cyclic loading must he able to account for changing boundary conditions 
during a specified load history. Usually, boundary conditions (free or fixed) 
in the finite— element method are satisfied by adding equations to, or deleting 
equations from, the overall system of equations. But the approach selected here 
was to connect two springs to each boundary node, as shown in Figure 3* One 
spring was used to satisfy boundary conditions in the x-direction, and the other 
to satisfy conditions in the y-direction. Therefore, all nodes in the system 
had two degrees of freedom. For free nodes, the spring stiffness, k gx or k gy , 
was set eq ual to zero. For fixed nodes, the spring stiffness was assigned an 
extremely large value (10 7 times the modulus of elasticity of the plate material). 
The spring stiffness was added to the diagonal coefficient in the conventional 
elastic stiffness matrix. The use of springs to satisfy boundary conditions 
was selected because an efficient technique to modify coefficients of the elastic 
stiffness matrix was incorporated into the nonlinear analysis program. This 
technique involved modifying the coefficients of the Cholesky decomposition [10] 
of the elastic stiffness matrix. The number of computer operations required to 
compute the modified Cholesky factors was of the order of n (n was the total 
number of degrees of freedom). In contrast, the number of operations required 
to obtain the original Cholesky factors was of the order of n 5 . A detailed 
discussion of the coefficient-modification technique may be found in [8] or [10]. 
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The coefficients of the elastic stiffness matrix are obtained from 


[K ] “ S./t B ] T CD dV m + [K s ] (5) 

m *!* 7 

where [D ] is the elasticity matrix and the diagonal matrix [K ] contains the 
e » 

elastic stiffness of the springs connected to the boundary nodes. 

The procedure for treating the nonlinear material behavior in the presence 
of changing boundary conditions remains unchanged from that previously presented 
for an elastic-plastic structure, except that all nodal displacements along the 
crack line were monitored to determine whether the nodes are to be released 
(crack extends), to open (crack opening), or to close (crack closure). 

To extend the crack, the crack-tip node was arbltarily chosen to be released 
at maximum load (the stiffness of the boundary spring was set equal to zero and 
the stiffness matrix was updated) and the crack tip was advanced to the next 
node. During the crack extension, the nodal force carried by the crack- tip node 
was released. To insure that the stresses and total strain increments in the 
adjacent elements satisfied the yield condition and the Prandtl-Reuss flow rule, 
the iterative procedure, previously discussed, was used to redistribute the 
force previously carried by the broken node, and to stabilize the plastic-load 
vector. 

During each load increment, the nodal displacements along the crack line 
were monitored to determine whether the crack surfaces had opened (positive 
displacement) or closed (negative displacement). If the crack surfaces had 
opened, the stiffness of the boundary spring was set equal to zero and the stiff- 
ness matrix was updated. If the crack surfaces had closed, the spring stiffness 
was set to the extremely large value and the stiffness matrix was updated. 



APPLICATION OF THE FINITE-ELEMENT ANALYSIS TO CYCLIC CRACK EXTENSION 

Elber [4], on the basis of fatigue crack-closure experiments with constant- 
amplitude loading, proposed the following equation for fatigue-crack-propagation 
rates 



(6) 


where C and n are material constants and ^ s the effective stress- 

intensity factor range- He proposed that the effective stress-intensity factor 
range be calculated by 


“eff “ V 55 a (7) 

where a is the half-length of the crack, a is a boundary-correction factor, 

and /CsS is the effective stress range, given by 
eix 




eff 


• S 


max 


S for S 
o o 


> S , 

— min 


( 8 ) 


where S is the maximum stress, S , is the minimum stress, and S is the 
max min ° 

crack-opening stress. Thus, the crack is assumed to propagate only during that 
portion of the load cycle in which the crack tip is open. Equations (6)- (8) 
suggest, then, that the crack opening stress can significantly influence crack 
growth under both constant- and variable- amplitude loading. 

The following sections give the results of applying the finite-element 
analysis to an extending crack under cyclic loading. The center-cracked panel. 
Figure 1, was subjected to either constant -amplitude or two-level block loading. 
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In this study, no attempt was made to establish a failure criterion for crack 
growth. For any cyclic loading, the crack- tip node (A, B, or F in Fig. 2) 

was arbitrarily chosen to be released at the maximum applied stress regardless 
of the magnitude of the applied stress and of any prior stress history. Thus, 
the model provides no direct information on the amount of crack growth per cycle; 
this information must be obtained from Equation (6). Instead, the analysis pro- 
vides only the crack opening stress, S q , to be used in Equation (8). 

The finite-element idealization and the coordinate system used for the 
center-cracked panel are shown in Figure 2. The initial crack tip was located 
at node A (see Fig. 2), so that the initial crack half-length, a^ was 2J.3 mm. 
The total panel width was h60 mm. The elastic-stress concentration (ratio of 
a in the most highly stressed crack- tip element to the applied stress) was 

yy 

approximately 7* 

The panel material was assumed to be elastic-perfectly plastic with a 

2 

tensile (and compressive) yield stress, a , of 350 MN/m and a modulus of elas- 
ticity of 70, 000 MN/m . These properties are representative of an aluminum alloy 
material. The cyclic stress-strain curve assumed for this material is shown in 
Figure 4. The cyclic stress-strain curve was also taken to be the stabilized 
curve; that is, the stress-strain behavior was assumed to be unaffected by fur- 
ther cycling. (Crews [11] has shown that in some materials the local stress- 
strain behavior at a notch root stabilizes in approximately 10 cycles. ) The 
solid line in Figure 4 shows the behavior under tensile loading. The dashed 
line indicates the typical behavior during unloading. The dash-dot line shows 
the behavior during reloading from a compressive plastic state. 
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Crack Extension Under Constant -Amplitude Loading 

Stress Level — The constant-amplitude loading applied to the center-crack 
panel is shown in Figure 5- The maximum gross stress was 0. 3 o q * As the cyclic 
stresses were applied, the crack initially opened at an infinitesimal stress due 
to the assumption of no prior plastic deformation. The symbol x indicates the 
point at which the most highly stressed crack-tip element initially yielded. At 
the maximum applied stress, node A (see Fig. 2) was allowed to displace (crack 
extension) and the crack tip advanced to node B. (Note that the crack extension 
increment, node A to node B, is solely a consequence of the element mesh size 
used and does not imply a crack-growth law. ) During unloading, node A was 
found to close at a positive stress (solid symbol) and was constrained against 
further motion. When the panel was reloaded, node A opened (open symbol) at a 
slightly lower stress than at which it had previously closed. Again, upon 
reaching the maximum stress, the next node (B) was also allowed to displace. 
During unloading, node (B) closed at a slightly higher stress than that at which 
node A had closed on the previous cycle. Further cycling and sequential release 
of nodes indicated approximately the same closure and opening stresses. The 
crack-opening stresses were approximately 30 percent of the maximum applied 
stress. 

Figure 6 shows another cyclic load history investigated, for which 

S = O.ko . Also shown are the calculated crack-closure and crack-opening 
max o 

stresses during constant-amplitude crack extension, indicated by the solid and 
open symbols, respectively. In this case, the contact (closure) stresses near 
node A after the first cycles were so large that material behind the crack tip 
yielded in compression and the subsequent opening stress was considerably lower 
than the previous closure stress. The closure stress rapidly stabilized to 
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approximately 1+0 percent of the. maximum applied stress. The opening stress 
increased after each cycle until it converged to the previously established 
closure stress. 

Because the crack-closure and crack-opening stresses are a function of the 

crack surface displacements, the variation of surface displacements with crack 

extension was investigated. . The crack surface displacements for constant- 

amplitude crack extension with S = 0. kcr are shown in Figure 7* The dis- 
r max o 

placement, v, in the y-direction (Fig. 2) is plotted as a function of the 
coor din ate location, x. The lowest curve shows the displacements at maximum 
stress with the crack tip located at node A (before crack extension). The other 
curves represent the crack-opening displacements at the maxi mum applied stress 
after each increment of crack extension. The sharp knee (slightly to the left 
of node A) in the displacement curves is an indication of the extent of severe 
plastic deformation. 

To gain a qualitative understanding of the relationship between the residual 

plastic deformations near the crack tip and the crack-closure stresses, crack- 

surface displacements were computed for two similar configurations. The first 

configuration is a crack, with the tip initially at node E, being opened for the 

first time by the maximum applied stress (S = 0. W ). The second configura- 

max o 

tion is a crack which was grown incrementally from node A to node E by cyclic 
loading (also with S = O.ko ). The results are shown in Figure 8. The 
dash-dot curve is the surface displacements for the first configuration under 
the maximum stress. The solid curve is for the second configuration under the 
same maximum stress. The f ully dashed curves show the displacements for the 
second configuration under a sequence of lower stresses for which the crack has 
closed successively from node E to node A. 
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The difference between the maxi mum surface displacements for the two config- 
urations, as identified by the shaded region in Figure 8, is a measure of the 
residual plastic deformations which are left in the wake of the advancing crack 
tip, and which Fiber has postulated as the major cause of crack closure during 
unloading. 

The displacement curves for lower stress levels further illustrate the 
crack-closure phenomenon* For example, a drop to 0.17^ Q applied stress 
caused the crack to close back to node D. Upon complete removal of the applied 
stress, the crack closed back to node A, but remained open to the left of node A. 

The stress distributions near the crack tip associated with these stages of 
unloading are shown in Figure 9 * At the maximum applied stress, the near-tip stress 
a (solid curve), reaches at plateau in front of the crack tip (characteristic of 

yy 

an elastic-perfectly plastic material). Because the local stresses are computed 
at the centroid of the elements (®= 0. 3mm from the crack surface), the stresses do 
not drop to zero immediately behind the crack tip. During unloading, the crack sur- 
faces contact initially at node D and the crack surfaces between nodes D and E begin 
to support compressive stresses, as illustrated by curves 3-5 of Figure 9 • 

As further illustration of the influence of stress level on crack closure. 
Figures 10 shows the results of constant -amplitude crack extension with 

S a o .5o. In contrast to the behavior at lower stress levels, this case 
max * o 

required six crack— growth increments to stabilize the crack— opening stress. 

Again, the contact stresses near nodes A, B, and C caused yielding in com- 
pression, as evidenced by the difference between the crack closure and opening 
stresses. In this case, the closure stress required a few cycles to stabilize. 
The opening stress increased with each cycle until it converged to the stabi- 
lized closure stress at approximately 50 percent of the maximum stress. 
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Thus, for the cases considered here the magnitude of the maximum applied 
stress apparently established the crack-closure stress which was found to be 
higher for higher values of maximum applied stress. In addition, with a suf- 
ficient number of cycles, the opening stress converged to the closure stress. 

Stress Ratio. - To study effects of stress ratio on crack closure, calcula- 
tions were made also for R = 0. 5 and R = -1. Figure 11 shows a cyclic load 

history with E s 0.5 and S =0.3 a • Again, nodes A, B, C, and D were 

max o 

released at succeeding occurrences of maximum stress* But in this case, the 
crack surfaces never closed* The solid symbols denote that the closure stress 
would have been approximately 0. 08q had the panel been allowed to unload 
beyond the minimum stress* If crack propagation rates are related to the effec- 
tive stress-intensity range (Eq. (6)), then, for stress ratios greater than 

approximately 0.3 and S = 0*3^ ? the present results suggest that crack- 

max o 

growth rates should be independent of the closure or opening stresses. This 

inference is consistent with the experimental results obtained by Katcher [12] 

on 2219-T851 aluminum alloy where the crack-closure stress was found to be lower 

than the minimum stress when the stress ratios exceeded approximately 0. 3. 

Figure 12 shows the crack- closure and crack-opening stresses for B = -1 

and S =0.3 0 * During the compressive loading, all nodes along the crack 

max o 

line closed. The crack-opening stress (at node A) after the first cycle was 
considerably lower than the opening stress obtained for the R = 0 case (Fig* 5)* 
This implies that the applied compressive stress caused the material near node A 
to yield in compression. However, after further crack extension, the closure 
and opening stresses stabilized at approximately the same level that had been 
observed previously at R = 0 (Fig. 5)* Therefore, according to Equation (6) 
the fatigue crack-growth rates for R * 0 and R - -1 should be approximately 



the same. This inference is consistent with results obtained by" Hudson [13], 

who found very little difference between the crack-growth rates at R = 0 and 

R = -1 in test of two different aluminum alloys. 

Crack Extension Under Two-Level Block Loading 

High-to-Low Loading — The high-to-low block loading is shown in Figure 13* 

The maximum stress S (first level) was O.l+cr and S (second level) was 

max, o max 


0-3cf • The crack-closure stresses for the three cycles of S 


were identical 


to those previously shown in Figure 6. At the maximum stress, S , node D 

dctXp 

was allowed to displace. During unloading, node D closed at a stress which was 
approximately 1+0 percent higher than the stabilized opening stress for S max 
(dashed line). When the panel was reloaded, node D opened at about the same 
stress at which it had previously closed. During further cycling, the closure 
and opening stresses decreased after each cycle. If more than three cycles of 


S were to be applied, the closure and opening stresses would be expected to 

max p 

converge to the dash-dot line (the stabilized opening stress for S ; see 
^ max^ 


Fig. 5) as the crack grew out of the material yielded by S 


Because the 


opening stresses for S are considerably higher than the stabilized open- 

m&Xp 

ing stress for S the crack-grovth rates (computed from Eq. (6)) are 

in&Xp 

expected to be considerably lower than if the three cycles s max had not been 
applied. This behavior has been experimentally observed and is referred to as 

v. 

crack-growth retardation or delay (see, for example [14] ). 

Low-to-High Loading — Figure 14 shows the lov-to-high block loading, for 


which S = 0.3a and S = O.ta . The crack closure and opening stresses 
max^ o maXg o 

for S were identical to those shown in Figure 5* When the stress first 


reached S , node D was released. During unloading, node D closed at approxi- 

maxg 

mately kO percent of the maximum applied stress. At zero applied stress, the 
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material near node D yielded in compression and the subsequent opening stress 

was lower than the previous closure stress. During further cycling, the closure 

stress remained nearly constant while the opening stress increased after each 

cycle. The opening stress had nearly converged to the stabilized opening stress 

after only three cycles of S , similar to the behavior observed in Figure 6. 

maxg 

Because the opening stresses for S sire lower than the stabilized opening 

max^ 

stress, the crack-growth rates (computed from Eq. (6)) are expected to be higher 

than the crack-growth rates computed using the stabilized opening stress for 

S . This behavior in crack-growth rates has also been experimentally observed 
max 2 

and is referred to as crack-growth acceleration [ 15 ]. 

CONCLUDING REMARKS 

A two-dimensional (plane stress ) finite element analysis which accounts for 
both nonlinear material behavior and changing boundary conditions under cyclic 
loading has been presented and used to study crack extension and crack-closure 
in a center-cracked panel. The finite-element model of the panel was assumed 
to be composed of an elastic-perfectly plastic material and was subjected to 
either constant-amplitude or two-level block loading. 

Although the failure criterion used in growing the crack was arbitrary, the 
calculated crack-closure and crack-opening stresses under both constant-amplitude 
and two-level block loading were qualitatively consistent with experimental 
observations. The calculated crack-opening stresses, when used with Elber f s 
crack growth equation, gave crack -growth-rate trends (retardation or acceleration) 
consistent with experimental measurements. Therefore, the finite-element analysis 
gives further insight into the mech ani sm of fatigue-crack growth during cyclic 
loading. 
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Figure 2. Finite-element idealization of the center-crack panel. 
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Figure 3. Springs connected to boundary nodes of a finite- 

element mesh. 
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Figure 4. Cyclic stress-strain curve for an elastic-perfectly 

plastic material. 
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Figure 5. Constant-amplitude crack extension with 









Figure 7. Crack surface displacements under constant-amplitude crack extension 

with S * 0.4 a and R = 0* 
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Figure 9. Local crack— tip stress distributions during unloading after constant- 

amplitude crack extension (S = 0.4 q ) with crack tip at node E. 
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Figure 10. Constant-amplitude crack extension with S max = 0.5 a Q and R = 0 
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Figure 13. High-to-low block-program loading crack extension with °o 

and S - 0.3 0. 
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